// simply supported beam by Beam2D element
/*
*                ^ force or moment
*        e1      |
*   |-------------
*   p1           p2
*/

#include <gtest/gtest.h>
#include <gmock/gmock.h>
#define private public
#define protected public
#include <element.h>
#undef private
#undef protected

using namespace std;

inline void printStdArray6d(const StdArray6d &arr)
{
    for (int i = 0; i < 6; i++)
    {
        cout << arr[i] << ", ";
    }
    cout << endl;
}

inline void printStdArray6d(const StdArray6d* arr)
{
    for (int i = 0; i < 6; i++)
    {
        cout << (*arr)[i] << ", ";
    }
    cout << endl;
}


int main()
{
    // create particles
    Particle* p1 = new Particle(1, 0, 0);
    Particle* p2 = new Particle(2, 10, 0);

    // create section
    Rectangle* sect = new Rectangle(1, 2, 4);

    // create material
    LinearElastic* mat = new LinearElastic(1);
    mat->setDensity(10);
    mat->setE(1000);
    mat->setNu(0.31);

    // create elements
    Beam2D* e1 = new Beam2D(1, p1, p2, mat, sect);

    // create DOF
    DOF d1 = {.key={true, true, false, false, false, true},
              .val={0, 0, 0, 0, 0, 0}};

    p1->constraintDof(d1);

    // force
    StdArray6d force = {0, 0, 0, 0, 0, 0};
    // force[5] = 10;
    force = {10, 10, 10, 10, 10, 10};

    // solving parameters
    double h = 0.001;
    double zeta = 0.5;
    int nStep = ceil(50 / h);

    for (int i = 0; i < nStep; i++)
    {
        if (i == 0)
        {
            p2->setForce(force);
            p1->solveFirstStep(h, zeta);
            p2->solveFirstStep(h, zeta);
            e1->setParticleForce();
        }
        else
        {
            p1->solve(h, zeta);
            p2->solve(h, zeta);
            p1->clearForce();
            p2->clearForce();
            p2->setForce(force);
            e1->setParticleForce();
        }
        p1->constraintDof(d1);
    }
    // results
    cout << "------- particle force ------" << endl;
    cout << "p1: ";
    printStdArray6d(p1->force());
    cout << "p2: ";
    printStdArray6d(p2->force());

    cout << "------- particle displace ------" << endl;
    cout << "p1: ";
    printStdArray6d(p1->displace());
    cout << "p2: ";
    printStdArray6d(p2->displace());


    cout << "------- element force ------" << endl;
    cout << "e1, fx, sfy, sfz, mx, my, mz: " << endl;
    cout << "emfi: " << e1->emfi_.fx << ", " << e1->emfi_.sfy << ", ";
    cout << e1->emfi_.sfz << ", " << e1->emfi_.mx << ", " << e1->emfi_.my;
    cout << ", " << e1->emfi_.mz << endl;
    cout << "emfj: " << e1->emfj_.fx << ", " << e1->emfj_.sfy << ", ";
    cout << e1->emfj_.sfz << ", " << e1->emfj_.mx << ", " << e1->emfj_.my;
    cout << ", " << e1->emfj_.mz << endl;

    cout << "e1: " << e1->force_.transpose() << endl;

    delete p1, p2;
    delete mat, sect;
    delete e1;
}

/* result summery
* force = {10, 0, 0, 0, 0, 0};
* e1: emfi = {-9.99996, 0, 0, 0, 0, 0},
*     emfj = { 9.99996, 0, 0, 0, 0, 0}
*
* force = {0, 10, 0, 0, 0, 0};
* e1: emfi = {-9.26087e-07, -9.99997, 0, 0, 0, -99.9548},
*     emfj = {9.26087e-07, 9.99997, 0, 0, 0, -5.72348e-06}
*
* force = {0, 0, 10, 0, 0, 0};
* e1: emfi = {0, 0, 0, 0, 0, 0},
*     emfj = {0, 0, 0, 0, 0, 0}
*
* force = {0, 0, 0, 10, 0, 0};
* e1: emfi = {0, 0, 0, 0, 0, 0},
*     emfj = {0, 0, 0, 0, 0, 0}
*
* force = {0, 0, 0, 0, 10, 0};
* e1: emfi = {0, 0, 0, 0, 0, 0},
*     emfj = {0, 0, 0, 0, 0, 0}
*
* force = {0, 0, 0, 0, 0, 10};
* e1: emfi = {-1.99354e-08, 4.12885e-06, 0, 0, 0, -9.99996},
*     emfj = {1.99354e-08, -4.12885e-06, 0, 0, 0, 10}
*
* force = {10, 10, 10, 10, 10, 10}; // large deformation, sect(2, 4)
* e1: emfi = {-9.99996, -9.99993, 0, 0, 0, -106.573},
*     emfj = {9.99996, 9.99993, 0, 0, 0, 9.99999}
*
* force = {10, 10, 10, 10, 10, 10}; // small deformation, sect(20, 40)
* e1: emfi = {-9.99996, -9.99998, 0, 0, 0, -110.001},
*     emfj = {9.99996, 9.99998, 0, 0, 0, 9.99987}
*/